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We propose and analyze an algorithm to approximate distribution func- 
tions and densities of perpetuities. Our algorithm refines an earlier approach 
based on iterating discretized versions of the fixed point equation that defines 
the perpetuity. We significantly reduce the complexity of the earlier algo- 
rithm. Also one particular perpetuity arising in the analysis of the selection 
algorithm Quickselect is studied in more detail. Our approach works well for 
distribution functions. For densities we have weaker error bounds although 
computer experiments indicate that densities can also be approximated well. 
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1 Introduction 

A perpetuity is a random variable X in R that satisfies the stochastic fixed-point equa- 
tion 

X = AX + b, (1) 

where the symbol = denotes that left and right hand side in ([!]) are identically distributed 
and where (A, b) is a vector of random variables being independent of X, whereas de- 
pendence between A and b is allowed. 
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Perpetuities arise in various different contexts: In discrete mathematics, perpetuities 
come up as the limit distributions of certain count statistics of decomposable combina- 
torial structures such as random permutations or random integers. In these areas, perpe- 
tuities often arise via relationships to the GEM and Poisson-Dirichlet distributions; see 



Arratia, Barbour, and Tavare (20031 for perpetuities, GEM and Poisson-Dirichlet dis- 
tribution in the context of combinatorial structures; see Donnelly and Grimmett (19931 
for occurrences in probabilistic number theory. In the probabilistic analysis of algo- 
rithms, perpetuities arise as limit distributions of certain cost measures of recursive 



algorithms such as the selection algorithm Quickselect, see e.g. Hwang and Tsai (20021 
or Mahmoud, Modarres, and Smythe (19951. In insurance and financial mathematics, 



a perpetuity represents the value of a commitment to make regular payments, where b 
represents the payment and A a discount factor both being subject to random fluctua- 



tion; see, e.g. |Goldie and Maller| ( |2000| ) or |Embrechts, Kliippelberg, and Mikosch| ( |1997 
Section 8.4). 

As perpetuities are given implicitly by their fixed-point characterization Q, properties 
of their distributions are not directly amenable. Nevertheless, various questions about 
perpetuities have already been settled. Necessary and sufficient conditions on (A, b) 
for the fixed-point equation ([!]) to uniquely determine a probability distribution are 
discussed in Vervaat ( 1979| ) and |Goldie and Mailer (20001. The types of distributions 
possible for perpetuities have been identified in Alsmeyer, Iksanov, and Rosier (2007 1 . 



Tail behavior of perpetuities has been studied for certain cases in Goldie and Griibel 



(1996). 



In the present article, we are interested in the central region of the distributions. The aim 
is to algorithmically approximate perpetuities, in particular their distribution functions 
and their Lebesgue densities (if they exist). 



For this, we apply and refine a method proposed in Devroye and Neininger (2002) that 



was originally designed for random variables X satisfying distributional fixed-point equa- 
tions of the form 

r=l 



(2) 



where . . . , X^ K \ (Ai, . . . , Ak, b) are independent with X^ being identically dis- 
tributed as X for r = 1, . . . , K and random coefficients A±, . . . , Ak, b, and K > 2. 

The case of perpetuities, i.e., K = 1, structurally differs from the cases K > 2: The 
presence of more than one independent copy of X on the right hand side in ([2]) often 
has a smoothing effect so that under mild additional assumptions on (A\, . . . , Ak, b) 
the existence of smooth Lebesgue densities of X follows, see Fill and Janson ] poooj ) 
and Devroye and Neininger (2002 1 . On the other hand, the case K = 1 often leads to 



distributions C(X) that have no smooth Lebesgue density; an example is discussed in 
Section [5] 
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Our basic approach to approximate perpetuities is as follows: A random variable X 
satisfies the distributional identity ([TJ if and only if its distribution is a fixed-point of 
the map T on the space A4 of probability distributions, given by 

T:M^M, /i ■-> C{AY + b), (3) 

where Y is independent of (A, b), and C(Y) = ^. Under the conditions \\A\\ p < 1 
and ||&||p < oo for some p > 1, which we assume throughout the paper, this map is a 
contraction on certain complete metric subspaces of M. Hence, C(X) can be obtained 
as limit of iterations of T, starting with some distribution jj,Q. 

However, it is not generally possible to algorithmically compute the iterations of T 
exactly. We therefore use discrete approximations (A^ n \b^) of (A, b), which become 
more accurate for increasing n, to approximate T by a mapping T^ n \ defined by 

f ( ft ) : M -» M, fJ, h-> C[A^Y + , 
where again Y is independent of (A^ n \ b^" 1 ') and C(Y) = fj>. 

To allow for an efficient computation of the approximation, we impose a further discreti- 
sation step (-) n , introduced in Section [2] defining 

: M -» M, n ^ ^(lA^Y + 6 (n) ) ) , 

where Y is independent of (A^ n \ b^) and C(Y) = /j,. 

In Section 2 we give conditions for o T^" 1 ) o • • • o T^\^ ) to converge to the 
perpetuity given as the solution of ([T]). To this aim, we derive a rate of convergence in 
the minimal L p metric £ p , defined on the space M p of probability measures on M. with 
finite absolute pth moment by 

£ p {u,^:=M{\\V-W\\ p : C(V) = u, C{W) = /i} , for^,/iG7W p , (4) 

where ||-|| denotes the L p -norm of random variables. To get an explicit error bound 
for the distribution function, we then convert this into a rate of convergence in the 
Kolmogorov metric g, defined by 

g{v,fi) := sup \F v {x) - F^{x) \ , 

where F u , denote the distribution functions of v, \i G Ai p . This implies explicit rates 
of convergence for distribution function and density, depending on the corresponding 
moduli of continuity of the fixed-point. 

For these moduli of continuity we find global bounds for perpetuities with 6 = 1 in 
Section |4} For cases with random b, we have to derive these moduli of continuity indi- 
vidually. One example, connected to the selection algorithm Quickselect, is worked out 
in detail in Section [U 
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We analyze the complexity of our approach in Section [3] As a measure for the com- 
plexity of the approximations for distribution function and density, we use the number 
of steps needed to obtain an approximation that has distance, in supremum norm, of at 



most 1/n to the true function. Although we generally follow the approach in Devroye 



and Neininger (2002), we can improve the complexity significantly by using different 



discretisations. For the approximation of the distribution function to an accuracy of 1/n 
in a typical case, we obtain a complexity of 0(n 1+£ ) for any e > 0. In comparison, the 
algorithm described in Devroye and Neininger ( |2002 ), which originally was designed for 



fixed-point equations of type ([2]) with K > 2, would lead to a complexity of 0{n i+£ h 
if applied to our cases. For the approximation of the density to an accuracy of we 
obtain a complexity of 0(n 1+1 / Q+e ) for any e > in the case of Q-H61der continuous 
densities, cf. Corollary |3.2| 



An extended abstract of this article appeared in Knape and Neininger (2007 



2 Discrete approximation and convergence 



Recall that our basic assumption in equation ([I]) is that ||-A|| p < 1 and ||6|| < oo for 
some p > 1. To obtain an algorithmically computable approximation of the solution of 
the fixed-point equation ([T|, we use an approximation of the sequence defined as follows: 
We replace (A,b) by a sequence of independent discrete approximations (A^ n \b^), 
converging to (^4, b) in pth mean for n — ► oo. To reduce the complexity, we introduce a 
further discretisation step (•} , which reduces the number of values attained by X n : 



X :=(EX) , I„:=i("'l„. 1 + & W , X n := (X n ) n , n > 1. 
We assume that the discretisations A^ n \ b^ and (■) satisfy 



- A 



b {n) - b 



<R b (n), (X n ) n -X n <R x (n 



(5) 



(6) 



for some error functions Ra, Rb and Rx, which we specify later. 
Furthermore, we assume that there exists some £ p < 1, such that for all n > 1, 

p 

which in applications is easy to obtain, since \\A\\ p < 1. 



(7) 



By arguments similar to those used in Fill and Janson (2002) and Devroye and Neininger 



(2002) we obtain the following convergence rates for the approximations X n to converge 



to the corresponding characteristics of the fixed-point X. We use the shorthand notation 
£ P (X,Y) :=£ p (£(X),C(Y)). 
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Lemma 2.1. Let (X n ) neNo be defined by §5§ and £ p as in ^fj. Then 

n-l 

£ P (x n ,x) <t; \\x-x Q \\ p + Y,SiRin-i), 

i=0 

where R(n) := Rx{n) + -Ra(^) ||X|| + Rb(n) for the error functions in 
Proof. We have 

U x n, X) < £ p (X n , X n ) + £ p (X n , X) 



< 



, x„ 



Xr 



+ e p (x n ,x). 



The first summand is bounded by ^ and for the second summand we have 



i P (X n ,X) < 



x n — X 



A^X n ^x + b {n) -AX-b 



< 



< 



A^X n ^ - AX 



+ 



b {n) - b 



A^ n \X n ^ - X) — (A — A^)X 



+ 



A {n) 



n-l — X\\ p + 



A - A^ 



b (n) - b 

x\\ + 



v 

b^ - b 



(8) 



(9) 



where in the last step we use that A^ and (X n -± — X) as well as (A — A^ n ') and X are 
independent by assumption. 

Now we use that the infimum in the definition of l p in Q is attained and assume 
additionally, that X n -\ and X are chosen with — X\ p = £ p (X n -i, X). Combining 

this with Q and using the bounds given in Q and 0, we obtain 



ip(X n ,X) < R x {n) + e p e p (X n ^,X) + R A {n) \\X\\ p + R h {n), 
and the claim then follows by induction. 



□ 



To make these estimates explicit we have to specify bounds for Ra(u), Rb(n), and 
Rx(n). We do so in two different ways, one representing a polynomial discretisation of 
the corresponding random variables and one representing an exponential discretisation. 
Better asymptotic results are obtained by the latter one. 

Corollary 2.2. Let X n ,n 6 No be defined by (|5| and £ p as in ([7]), and assume 

R A {n) < C A 4> Mn) < C b — r , R x {n) < C x ^, 
n r n r n r 

for some r > 1. Then, we have 



?p{X n , X) < c r — , 
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where 



C r 



r r \\X-X n 



e log (1/6 



+ 



(c x + c b + c A \\X\ 



r+1 



(10) 



Proof. Using Lemma 2.1 we get 



n-l 



,(X n ,X) < % \\X - X \\ p + {C x + C A \\X\\ p + C b ) £ 



i=0 



(n — i) 1 



(11) 



For the first summand, we use that the function x i— > x r £p has its maximum at x = 

?7iog(i/g. 

To see that the second summand is of order n~ r , note that l/(n — i) < (i + l)/n for all 
n > 1 and < i < re — 1. This implies that for £ p < 1, 



n— 1 1 n— 1 

y^W<-y(^ + i) r e 



i=0 



(re — i) r n 



i=0 



i=0 

7*' 



(i-g' 

where the last equality is obtained by differentiating the geometric series r times. □ 



Remark 2.3. In Corollary 2.2 we are merely interested in the order of magnitude of 
£ p (X n ,X) without a sharp estimate of the constant C r . When evaluating the error in 



an explicit example, we can evaluate (11) directly to obtain sharper estimates 



Corollary 2.4. Let X n ,n 6 No be defined by ^ and £ p as in (J7|, and assume 

R A (n)<C A ± R b ( n )<C b ± Rx(n)<C x ± 
for some 1 < 7 < l/£„. Then, we have 



£ P (x n , X) < C 7 — , 



where 



C 7 :— \\X — Xo\\ p + 



(C x + C b + C A \\X\\ p 



(12) 
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Proof. Using Lemma 2.1 we get 



n-l 



i P (x n ,x) < e; \\x - x \\ p + (c x + c A \\x\\ p + c fe ) 7 - n ^e;y 

i=0 



(13) 



and the assumption on 7 implies that both summands are 0(7 n ) with the constant 
given in the lemma. □ 



Then, the distance in the Kolmogorov metric can be bounded by 



Lemma 2.5. Let X n and C r be as in Corollary \2.2\ and X have a bounded density fx- 

(14) 



e(X n ,x)< (O-b+i^ll/xlLj 



\p/(p+i) 



n 



-rp/(p+l) 



Similarly, for X n and C 7 as in Corollary 2.4, we have 



g(X n ,X) < (c r ( P+ l)^\\f x \Q P/[P+l) ^^. 



(15) 



Proof. We use Lemma 5.1 in Fill and Janson (20021, which states, that for X with 
bounded density fx and any Y, 

q(Y,X)< ((p + l) 1/p \\fx\L £ P (Y,X)Y /iP+1) forp>l. 



Using Corollaries 2.2 and |2.4| respectively, we get the stated result. 



□ 



Remark 2.6. In some cases, we can give a similar bound, although the density of X is 
not bounded or no explicit bound is known. Instead, it is sufficient to have a bound for 



the modulus of continuity of the distribution function Fx of X, cf. Knape (2006). 



To approximate the density of the fixed-point, we define 

F n (x + 5 n ) - F n (x - S n ) 



fn{x) 



25 n 



(16) 



where F n is the distribution function of X n . For this approximation we can give a rate 
of convergence, depending on the modulus of continuity of the density of the fixed-point, 
which is defined by 



A fx (5):= sup \fx(u)-fx(v)\, 5>0. 

\u—v\<8 



Lemma 2.7. Let X have a density fx and let X n ,n € No be defined by Then, for 
f n defined by (16) and all 5 n > 0, 



\fn-fx\L < T g(X n ,X) + A fx (5 ri 

On 
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Proof. For any x, we have 



\fn(x) ~ fx(x)\ 



< 



F n (x + 6 n ) - F n (x - 5 n ) F{x + 8 n ) - F(x - 5 n ) 



+ 



25 n 

F(x + S n ) - F(x - 6 n ) 



2S n 

< ~ g(X n ,X) + \fx(x + y)~ /; 

On 26 n J_ Sn 

< ^g(X n ,X) + ±- [ \ x {y)dy- 

On n Jq 

The assertion follows since Af x is monotonically increasing. 



25 n 
fx(x) 



+ 



(x)\ dy 



□ 



Corollary 2.8. Let X have a bounded density fx, which is Holder continuous with 



exponent a £ (0,1]. For polynomial discretisation X n and C r as in Corollary 2.2 and 

«5„:=Ln-^/« Q+1 )( p+1 )) 



f n defined by (16) with 



with an L > 0, we have 



\p/(p+i), 



n-fx\L< [{Cr(p+l) 1/P \\fx\U ' '/L + CL«\ n -« r P/(( Q + 1 )(P+l))_ 



For exponential discretisation X n and C 7 as in Corollary 2.4 and f n defined by (16) 
with 

5n . = ^ 7 -pn/((a+l)(p+l)) ) 



with an L > 0, we obtain 

11/n-ZxlL < ((C 7 (p+l) 1 / p ||/x|| 00 ) P/(P+ % + cL^ 7 W((a+Db+D). 



Remark 2.9. If X is bounded and bounds for the density fx and its modulus of 
continuity are known explicitly, the last result is strong enough to construct a perfect 
simulation algorithm based on von Neumann's rejection method. Corollary |2.8| can 



be turned into such an algorithm as done in Devroye (2001) for the case of infinitely 



divisible perpetuities with approximation of densities by Fourier inversion, |Devroye, Fill, 
and Neininger ( 2000 ) for the case of the Quicksort limit distribution and Devroye and 



Neininger| ( |2002 ) for more general fixed-point equations of type ^ . 



3 Algorithm and Complexity 

In this section, we will give an algorithm for an approximation satisfying the assumptions 
in the last section for many important cases. We assume that the distributions of A and 
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b are given by Skorohod representations, i.e. by measurable functions <p, tp : [0, 1] — * M, 
such that 

A = ip(U) and b = ip(U), (17) 

U being uniformly distributed on [0,1]. Furthermore, we assume that IMloo < 1 an d 
that both functions are Lipschitz continuous and can be evaluated in constant time. 
Now we define the discretisation (•) by 

(Y) n := [s(n)Y\/s(n), (18) 

where s(n) can be either polynomial, i.e. s(n) = n r or exponential, s(n) = r y n . Defining 

A^:=v{{U) n ) and 
b (n) :=V)(( m ), 



the conditions on ip and ip ensure that Corollary |2.2| and 2.4 can be applied. 
We keep the distribution of X n in an array A n , where 

An[k] := F[X n = k/s(n)\ 

for k £ Z. Note however, that as A and b are bounded, »4n[A;] = at least for \k\ > 
s(n)Q n , where Q n can be computed recursively as Q n = [H^H^ Qn-i + ll^llool an< ^ 
Qo=\\\X \U=\EX]. 

For simplicity we assume that s(0) = s(l) = 1 and that s(n) E N for all n. The core of 
the implementation is the following update procedure: 

procedure update(A„_i, A n ) 
for i <— to s(n) — 1 do 

for j < s(n — 1) Q n -i to s(n — 1) Q n -i do 



s(n) 



s(n) y?(ti) 



s(n — 1) 



end for 
end for 
end procedure 

Furthermore, we use a procedure INITIALIZE(^4 n , ri), which creates A n as vector with 
2s(n)Q n components with A n [k] = for —s(n)Q n < k < s{n)Q n . 

The whole algorithm then looks like this: 
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initialize(^4 05 0) 



s(0) EX 



1 



(19) 



for n <— 1 to N do 

INITIALIZE(Ai, n) 
UPDATE (An-i,An 

end for 

return An 



Note, that (19) determines that we start the approximation with Xq as defined in 

The complete code for polynomial discretisation for the example in Section [5] imple- 



mented in C++, can be found in Knape (2006). 



To approximate the density as in (16) with 5x = d/s(N) for some d G N, we compute a 
new array T>n by setting 

s{N) 



k+d 



V 



2d 



j=k-d+i 



To measure the complexity of our algorithm, we estimate the number of steps needed to 
approximate the distribution function and the density up to an accuracy of 1 /n. For the 
case that X has a bounded density fx which is Holder continuous, we give asymptotic 
bounds for polynomial as well as for exponential discretisation. We assume the general 



condition (17). 



Lemma 3.1. Assume that X has a bounded density fx, which is Holder continuous with 



exponent a G (0,1]. Using polynomial discretisation with exponent r, cf. Corollary 2.2 



we can calculate for any n G N approximations F, f of the distribution function F and 
the density f of X with 



1 



< 

oo n ' 



F-F 

in time Tp(n) and Tf(n) respectively with 



f-f 



1 

< - 
oo n 



Tp(n) = 0( n ( 2 + 2 /^+ 1 )/p) and T f (n) = o^i+i/^+Db+D/M) . 

Using exponential discretisation with parameter^ as in Corollary \2.J\ approximation to 
the same accuracy takes time 

T' F (n) = O (n^ p+1 ^ p log nj and T' f {n) = O (n^ 1+1 ^ p+1 ^ p log n) 

for the distribution function and the density of X respectively. 
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Proof. In one execution of \JPDATE(Ak-i, Ak) , the outer loop is executed s(k) times. 
The assumptions on A and b ensure that Qk = 0(k), so we have 0(ks(k)) runs of the 
inner loop and the whole procedure takes time 0(k s(/c) 2 ). Hence, for any N £ N, 
finding An costs time 

O^Y^ksik) 2 ^ =0(^N 2 s{N) 2 y (20) 
For discretisations with s(n) = n r we get a running time of 0(N 2r+2 ) to find An, and 



(14) in Lemma 2.5 ensures that for the corresponding distribution function Fn of Xn, 

WFn-FW^ < CN~ rp ^ p+1 \ 

Setting N = (Cn) {p+1)/{rp) and F := F N , we get an approximation of the stated accuracy 
in time 

T F (n) = 0(N 2r + 2 ) = o(n( 2+2 A'XP +1 )/P). 

For the density of X we use Corollary Q and N' = (c' n ) {a+1)ip+1)/{arp) to obtain the 
stated bound. 

When using exponential discretisation, s(n) = j n , we need time 0(N 2 j N ) to find An- 
Using the corresponding results in Lemma |2.5| and Corollary |2.8| ensures the stated 
running times. □ 



Corollary 3.2. Assume (|l7j) and that X has a bounded density fx, which is Holder 
continuous with exponent a £ (0, 1]. Then, using exponential discretisation as in Corol- 



lary 2.4. approximation to an accuracy of 1/n takes time 0(n 1+e ) for the distribution 



function and time 0{n l+l / a+£ ) for the density of X for all e > 0. 

Proof. Note that ||<£>|loo — 1 an d ||^4|| p < 1 for some p > 1 implies that ||^4|| p < 1 for all 



p > 1. Thus, in Lemma 3.1, p can be chosen arbitrarily large. □ 



4 A simple class of perpetuities 

In order to make the bounds of Section [2] explicit in applications, we need to bound 
the absolute value and modulus of continuity of the density of the fixed-point. For a 
simple class of fixed-point equations, we give universal bounds in this section. For more 
complicated cases, bounds have to be derived individually, which we work out for one 
example in Section [5] 

For fixed-point equations of the form 

X = AX + 1 with A>0, (21) 

where A and X are independent, we can bound the density and modulus of continuity 
of X using the corresponding values of A. 
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Lemma 4.1. Let X satisfy fixed-point equation (21) and A have a density /a- Then X 
has a density fx satisfying 



fxiu) 

and fx(u) = otherwise. 



*> 1 fu—l^ 

fA[ — )fx(x)dx, foru>l, 



l x 



x 



(22) 



Proof. From the fixed-point equation we can see that X > 1 almost surely. Now let Fx 
be the distribution of X. Conditioning on X, we get for any Borel set B: 

/oo 
F[Ax + 1 G B]dF x (x) 

fxA+\{u)dudF x {x) 
r u-\ 



1 JB 



i Jb x 



fA 



X 

u-1 



dudF x {x) 
T x (x) du, 



where we can use Fubini's theorem in the last step, because the integrand is product 
measurable. The claim follows, as this is just the definition of a Lebesgue density. □ 

Corollary 4.2. Let A have a bounded density fA- Then X has a density fx satisfying 

ll/xlloo < II/aIL • 



Proof. Using Lemma 4.1 we get 



11/xlL < II// 



E 



1 

X 



but X > 1 implies E[l/X] < 1, so the claim follows. 



□ 



Corollary 4.3. Let A have a density fA, and Af A be its modulus of continuity. Then 
the modulus of continuity A t x of fx satisfies 

A fx (S)<A fA (S), 5>0. 



Proof. Using (22), we obtain for any u, v G 



\f.x(tn- fx(r)\ < I -fx(x) 

ll X 



fA 



u-1 



fA 



v - 1 



dx. 



(23) 
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But x > 1 and the modulus of continuity Af A is monotonically increasing by definition, 
so we can bound 



fA 



u-1 



fA 



v - 1 



X 



< A 



fA 



u — V 



<A fA (\u-v\), 



and plugging this into inequality (23), we obtain 



\fx(u)-fx(v)\ <E 



1 

X 



A fA (\u-v\ 



Now we use that E[l/X] < 1 and take the supremum over all suitable u, v. 



□ 



This result is only useful if the density of A is continuous, but we can extend it to many 
practical examples, where /a has jumps at points in a set Ia- We use the jump function 
of /a, defined by 

J f A ( s ) = /a( s ) - lip/A^), s > 

xjs 

and a modification of Ja where we remove all jumps, 

fA ■= fA ~ ^ ^(^[s.oo)- 

seJ 4 \{0} 

Since X > 1, we now denote by Af x the modulus of continuity of the restriction of fx 
to (1, oo). 

Lemma 4.4. Let A have a bounded cddldg density fA- Then, for all 5 > 0, 

\JfA(*)\8 



±fM < ±hW + \\fx\L E 



seJ A \{0} 



Proof. We give the proof for the case that fA has only one jump, say in so > 0. The 
general case then follows similarly. For 1 < u < v, we have 



C°° 1 

\fx{u)-fx{v)\ < / -fx{x) 



fA 



u-1 



fA 



V - 1 



dx. 



We define 



u — 1 . , v — 1 
a:= VI, (5:= VI 



so so 

and divide the range of integration into the three intervals (1, a], [a, (3], and [/3, oo). Now, 
in the first and third interval, differences of values of fA and fA coincide. Moreover, for 
x £ [a, (3] we have 



fA 



u-i 



fA 



v-i 



< 



Ia 



u-1 



X 



Fa 



v - 1 



+ \ J f A ( s o)\ 
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Putting everything together we obtain 



\fx(v) - fx(u)\ < 



f°° 1 

< / -fx(x) 

rco -i 

< I ~fx{x) 

X 



Fa 
h 



u - 1 

X 

II - 1 

X 



h 
h 



v - 1 

X 

r - 1 

x 



dx+ -fx{x)\Jf A {s )\dx 



dx+\\f x 



ixj —T~ I J /a( s o)| • 



■so 



We now bound the latter integral by Aj A (v — u) as in Corollary 4.3 and the claim follows 
by taking the supremum over all v — u < 5. □ 



5 Example: Number of key exchanges in Quickselect 

In this section, we apply our algorithm to the fixed-point equation 

X = UX + U(1-U), 



(24) 



where U and X are independent and U is uniformly distributed on [0, 1]. This equation 
appears in the analysis of the selection algorithm Quickselect. The asymptotic distribu- 
tion of the number of key exchanges executed by Quickselect when acting on a random 
equiprobable permutation of length n and selecting an element of rank k = o{n) can be 



characterized by the above fixed-point equation, see Hwang and Tsai (2002). 



We use our algorithm to get a discrete approximation of the fixed point. The plot of 
a histogram, generated with 80 iterations of the algorithms using for the discretisation 
s(n) = n 3 , can be found in Figure [!} 

In the following, we specify how the bounds in Section [2] can be made explicit for this 
example. 



Lemma 5.1. Let X be a solution of (24). Then, we have < X < 1 almost surely, and 
the moments are recursively given by E[X°] = 1 and 



fc-i 



E 



X h 



^ j\{2k-j + 1)! 



k > 1, 



in particular, K[X] = 1/3. 



Proof. Both claims follow directly from the fixed-point equation in ( 24 ) , using that the 
solution is unique. To compute the moments, note that E[C/ fc (l — U) k ~^ is equal to the 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Figure 1: Histogram of approximation for X = UX + U(l — U) 

Beta function B(/c + 1, A; — j + 1), so we have 

k-l 



E 



X k 



k + 1 



k-l 



fc! k\(k-j)\ 



/>■ ^iK^-iO (2fc-j + i)! 



and the assertion follows. 



Lemma 5.2. Let X be a solution of (24). Then, for all k £ N and e > 



P[X > 1 - e] < 2 ( k2 - k )/ 4 £ k / 2 . 

Proof. Using that X is supported by [0, 1], it is easy to show that for all e 

P[X > 1 - e] = P[UX + U(l-U)>l-e] 

< P[X >l-2e]¥[U>l-y/e\, 

and this inequality can be translated into 

nx>l-2e]> P[X ^- £] . 



Applying (25) k times, we get 



1 > F[X > 1 - 2 K e] > 



P[X > 1-e] 

2k(k-1)/4 £ k/2' 



This implies the assertion. 
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Lemma 5.3. Let X be a solution of (24). Then X has a Lebesgue density f satisfying 
f(t) = Ofort<0 or t > 1 and 

f(t) = 2 [ g(x, t)f(x)dx + [ g(x, t)f(x)dx for t G [0, 1] , (26) 
Jp t Jt 

Pt:=2\ft — 1, q(x,t):= — , 



where 



Proof. Let Px be the distribution of X. Then we get for any Borel set B by conditioning 
on X as in the proof of Lemma |4.1| 



"[X G B] = P [UX + U(l -U)eB] 



[ ¥[Ux + U(l -U)eB] d¥ x {x) 
Jo 



JB 



<p x (t)dt d¥ x (x) 




B JO 



if x (t)d¥ x {x) dt 



where ip x is a Lebesgue density of (1 + x)U — U 2 . The last step is valid by Fubini's 



theorem as (x,t) h-> <p x (t) is product measurable, cf. (28). 
Hence, X has a Lebesgue-density f(x) satisfying 

f(t) = I <p x (t)f(x)dx. 
Jo 

To find (p x , we observe that (1 + x)U — U 2 < (1 + x) 2 /4 and get 
P[(l +x)U- U 2 <t) = 



(27) 



TT ^l+x-J{l+x) 2 -U TT ^ l + a; + ^(l + x) 2 -4t 
U < 1 or U > 



for t < 0, 

1+x- +x) 2 - At 



for < t < x, 



1 - + x) 2 -At for x < t < (1 + x) 2 /4, 

1 otherwise. 
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To get a density, we differentiate with respect to t and rewrite as a function of x yielding 

2 



<Px(t) = < 



^/{l + x) 2 -At 
1 



for 2 y/t -l<x<t, 
for t < x < 1, 



^(l + x) 2 -4t 
otherwise. 



Plugging this into (27) we get the stated integral equation. 



(28) 



□ 



Remark 5.4. The integral of g(x,t) with respect to x can explicitly be evaluated: 

g(x, t) dx = log (l + x + x/(l + x) 2 -4t) . (29) 



Remark 5.5. We will see in Lemma 5.7 that f(x) has a version that is continuous on 
[0, 1]. For this version we have 



/(0) = E 



1 



1 + X 



0.759947956... 



Proof. Using integral equation ( 26 ) we have 

/(0) 



1 



io l + x 

and by expanding the geometric series we obtain 

1 



f(x)dx, 



E 



l + X 



£(-l) fc E 



fc=0 



which we can calculate to any accuracy using for the fcth moments the formula given in 
Lemma 15.11 □ 



In order to use Lemma |2.5| to bound the deviation of our approximation, we need an 
explicit bound for the density of X. We derive a rather rough bound here and see later, 
that we can use the resulting bound from our approximation to improve it. 

Lemma 5.6. Let f be the density of X as in Lemma \5.3[ Then 

< 18. 



Proof. To get an explicit bound for t £ [0, 1] we simplify the integral equation and obtain 



f(t)<2 [ g(x,t)f(x)dx. 
Jpt 



(30) 
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We know f(t) for t < 0, and we can bound g(x, t), if x is bounded away from pt. Hence 
we split the integral into a left part for which we already have a bound for / and a right 
part, in which we can bound g. For any 7 G (p t , 1], we have 



f(t)<2 r g {x,t)dx + 2 [ g{x,t)f{x)dx, 
Jpt J 7 



(31) 



where in the second integral, we can use that g is decreasing in x for any fixed t and 
bound g{x,t) < 5(7, i). 

For t < 1/4, we can use that pt is negative, and set 7 = 0. So the first integral vanishes 
and only the second remains and we obtain 



f(t)<2 [ g(x,t)f(x)dx<2g(0,t) [ /( 
Jo Jo 



x)dx 



To go on, we set 7 = 7*:= {p t + t)/2 and get with (31 ) 

f{t) < 2 ^ fg{x, t)dx + 2 <7( 7t , t) I f(x)dx, 



pt 



it 



where \i t := sup{/(t) : r <E (pt,7t)}- 



We can calculate the first integral using the integral of g given in ( 29 ) , 

nt 

/ 9{x,t)d. 

J Vt 



ix = log 1 + 



(1 - Vi) 2 + (1 - Vi)V 1 + Gy/i + t 



and for the second integral, we obtain 



1 ,1 
f(x)dx < / f{x)dx 
"it Jpt 



Putting everything together we get 



4v^ 



X > 1 -2(1 -Vt 



f(t)<2fi t h(t) + 4 



'[X> 1-2(1- y/i)] 



(1 - \/i)\Jl + 6^ + * 



(32) 



--:h{t), (33) 



(34) 



For t = 1/4 we have 7^4 = 1/8, and /x 1/ / 4 < 2y2 by (32), so 



/(l/4)<4v / 2 log 1 + l^y 17 +-^=<7. 



(35) 
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From the integral equation we get for 0<s<i<l/4 

/(*) ~ /(«) = / {9(x, t) - g(x, s))f(x)dx+ 
Jo 

rt 

(g(x,t) - g(x,s))f(x)dx+ / g(x,t)f(x)da 

J s 



+ 
> o, 

so / is strictly increasing on [0,1/4]. Therefore, the bound for t = 1/4 extends to aU 
t G [0, 1/4] =: Iq. To go on, we recursively define bo := and 

2 

, »>1, 



and 



hk-i '■= [ h, 



h := 
h + bk+i 



1 + 



1 , 

-j Ofc+1 



, fe > 1. 



For each interval I n we find a corresponding bound M n for /, using that = bi-\ and 
therefore (pt,Jt) C I„_i U 7 n _ 2 for t E I n . 



Furthermore we get for 1 /4 < t < 1 by differentiating the function h defined in ( 33 ) 
h'(t) = c t 



d (1-Vi) 2 d (l-y^)yi + 6>A+t 
(it 4^/t dt 



where q > 1. But the first summand is negative and for the second observe that 



d_ 
dt 



(1 - V~t)yi + 6Vt + t 



[1 -Vi)(3 + Vi)-(l + 6Vi + t) 



2 y/iy/l + 6y/t + t 
\-kJi-t 



iVl + 6-s/i + t 



< 0, 



hence h{t) is decreasing. 



The second summand in ( 34 ) can be bounded using Lemma |5.2| with k = 2 yielding 
4 



>[X> 1-2(1- yft] <i P[X>l-2(l-y / t)] ^ 



2 (1 - V*) 



(36) 



(37) 



(1 - Vt)Vl + 6^ + * 
So for t £ I n = (q„, /3 n ] we have 

/(*) < M n := ^(a,,) max{M n _i,M n _ 2 } + 4V2 
Evaluating this we obtain 

M = 7, Mi = 13, M 2 = 17, M 3 = 18, M 4 = 17. 
But for t > 63 we have /i(i) < 2/7 so the sequence (M n ) n>0 is decreasing for n > 4. □ 
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Lemma 5.7. Let f be the density of X as in Lemma 5.3. Then f is Holder continuous 
on [0,1] with Holder exponent 1/2: 



!/(*)-/(*) I < 



vr 



for < s < t < 1. 



(38) 



Proof. Using the integral equation given in Lemma 5.3 we have 



\f(t)-m\ <2 



pi 



g(x,t)f(x)dx- / g(x,s)f(x)dx 



+ 



+ 



g(x,t)f(x)- / g(x,s)f(x)dx 



(39) 



With explicit calculations we find 



pi 



g(x,t)f(x)dx- / g(x,s)f(x)dx 



and 



g(x,t)f(x)dx- / g(x,s)f(x)dx 



< 4 



< 



For details see Knape (2006). 



□ 



Remark 5.8. The latter lemma cannot be substantially improved, as in t = 1/4, the 
density f{t) is not Holder continuous with Holder exponent 1/2 + e for any e > 0, see 



Knape (2006). 



6 Explicit error bounds for X = UX + U(l - U) 

We can now combine the bounds for the density and its modulus of continuity with 



Lemma 2.5 and Lemma 2.7 to bound the deviation of an approximation from the solution 
of the fixed-point equation. 



To approximate the density / we set 

7(o) 

fn(x) : 



F n (x + 5 n ) - F n (x - 5 n ) 



2<L 







for < x < S n , 

for 5 n < x < 1, 
otherwise, 



where /(0) is given in Remark 5.5 and F n denotes the distribution function of X n . 

For the values used for the plot in Figure [TJ i.e. s(n) = n 3 and iV = 80, we can apply 
Corollary |2.2| and obtain: 
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Corollary 6.1. We have g(X 80 ,X) < 1.162 • 1(T 4 , and ||/ 80 - /H^ < 0.931. Further- 
more, we can improve the bound of Lemma 5.6 and bound ||/|| < 3.561. 



Proof. We have Ca = Cfc = Cx = 1, hence combining Lemma 5.6 and Lemma 2.5 
obtain 



we 



// n-1 ci v \1>. 

Q(x n , x) < U% \\x\\ p + (2 + £ j (p + 1)^ ii/il . 



p/(p+l) 



i=0 v " '' ' J 

The moments of X can be computed using Lemma 5.1 and we set [U] n : = |_n 3 C/J /n 3 
hence 



? r = li^ P = (^) 1P . 

Optimizing over p for n = 80, r = 3, and || H— < 18 yields 

q{X S0 ,X) < 5.1842 • 10~ 4 (40) 

for p = 12. 



Using for /(0) the value given in Remark 5.5 we obtain for the density 

\\fn-f\L< ±g(X n ,X) + 9\\f\\ 00 yfa, 

On. 



and optimizing over S n , using for the Kolmogorov metric the bound in (40), yields 

ll/so -/IL< 4.512 
for <5go = 3.44 • 10~ 4 (averaging 352 values). 

We can now use this to improve our bound for H/H^: Reading off the maximal value of 
our approximation (H/solloo — 2.630), we can now bound 

||/||oc<ll/80|loc + 11/80 -/IL< 7.142, 

and this in turn enables us to improve our bounds for the approximation, leading to 
g{X 80 ,X) < 2.2085 • 10" 4 and ||/ 80 - /IL < 1-8331 for 5 80 = 3.6 • 10" 4 . Repeating this 
strategy a few times, we get the stated values for p = 13 and <5so = 3.7 • 10~ 4 (averaging 
378 values). □ 

Remark 6.2. Using the realistic (but yet unproven) bound of \\fl\~ < 2.7 would give 
g(X 80 ,X) < 8.9809- 10" 5 (p = 13) and ||/so - /IL < 0.7101. Hence, our approach 
works well for the distribution function. However, we cannot show strong error bounds 
for the approximation of densities with our arguments. 

However, in the next section we see that for another example the algorithm approximates 
the densities much better than the error bounds indicate. 

In Table [T] the resulting error bounds for several possible discretisations with similar 
running time can be found. 
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Discret. 


N 


q(X n ,X) 


opt. p 


s(N) 


n 


22000 


0.00178 


14 


22000 


n 2 


430 


0.00025 


16 


184900 


n 3 


80 


0.00012 


13 


512000 


4 

n 


30 


0.00050 


3 


810000 


1.5™ 


35 


0.00070 


3 


1456110 


1.7 n 


27 


0.00187 


2 


1667712 



Table 1: Bounds for g(X n ,X) for comparable total running times (about 20h on a lap- 



top computer each). The discretisations are according to Corollaries 2.2 and 



2.4 By s(N) the number of atoms of the discrete approximation is denoted, 



cf. Section [3l 

7 An experimental view on error bounds 



We now apply our algorithm to another fixed-point equation for which the solution is 
explicitly known. We can then compare the approximation of our algorithm with the 
true density and distribution function and evaluate the actual error to get an idea of 
the quality of the error bounds proven in Section [2j Further examples can be found in 
Knape ( 2006 ) . It appears that the error bounds in Section [2] are rather loose and that 
the approximation is much better than indicated by our bounds. 

In the analysis of certain random interval splitting procedures the following fixed-point 
equation characterizes the distribution of a point to which a random sequence of nested 
intervals shrinks: 

d 1 + U „ „!-£/ 



X 



X + G- 



2 2 ' 

where G, U, and X are independent, G is Bernoulli (1/2) distributed and U is uniformly 



distributed on [0,1], see 


Chen, Goodman, and Zame 


( 


.984), Chen, Lin, and Zame ( 


1981 


I, 


Devroye, Letac, and Seshadri 


( 


1986 


), and 


Neiningei 


' ( 


2001 


) for details of the interval 



splitting context. 

To approximate the fixed-point, we use a symmetric discretisation for (A, b) instead of 



(18), setting 



(U) n := (2 [s(n)U\+l)/2s(n) 



(41) 



and s(n) = n 3 . 

To compute the bounds as given in Section [2j we can set Ca = Cb = 1/4, £ p = \\A\ 
and A is uniformly distributed on [1/2, 1], so 



9P+1 _ i 

I 4 iip = _ 

1 llp 2P(p+l) 



for p£N. 
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It is known that X is beta(2, 2) distributed, so we have the moments: 



p-i 



Furthermore, X has the density fix) = 6x(l — x), so 
Lemma 2.5 and Corollary 2.2 to obtain 



1.5. We can now use 



q(X n ,X) < [ 1.5 (p + l) 1/p [\\A\\^\\X\\ p + 



5+11X1 



N-l 

-£ 

i=0 



lAII 



(N-iy 



p 

p+i 



For iV = 50 we minimize over p and get p m i n = 5 and 

g(X 50 ,X) < 0.001043. 



(42) 



As we know the limit distribution, we can read off the true error from the output of our 
simulation and find 

g{X 50 ,X) « 0.000012. 

It is quite exactly of the order expected for a discretisation of step size 1/n 3 . Note that 
when approximating a differentiable function by a step function, step size and derivative 
impose an unavoidable error. Comparing our approximation to a direct discretisation 
by a step function of the same step size, the deviation is at most 1.5 • 10~ 8 . 

Now we look at the density. The modulus of continuity of the density of the beta(2,2) 
distribution can be bounded by A/(e) < 6e for all positive e. So for the function fjy, 



which we get by averaging over 25 n as in (16 1, we get with Lemma 2.7 



||/v-/|L < j^g(X N ,X)+68 N . 



We evaluate for N = 50, use the bound in (42), and minimizing over 5^q we obtain 



ll/so - /|L < 0.1583 

for = 0.01318, so we take the average over 3 296 values. 
Reading off the true error from the simulation we obtain 

||(/n-/)l[0iH5;0.985]|L~ 0.0003 

and \ f n {x) — f{x)\ < 0.02 for x < 0.015 or x > 0.985. The larger errors at the boundary 
are caused by the averaging procedure used to obtain f n . 

Acknowledgements: We thank the referee for careful reading, pointing out some 
inacurracies and helping improve the presentation of the paper. 
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